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ABSTRACT 
Linear  state  signal  feedback  is  used  to  obtain  expomiential 
response  from  fourth  order  systems.   Characteristic  equation  roots 
are  selected  to  provide  the  desired  exponential  response  with  con- 
straint on  initial  conditions  and  system  acceleration „  A  digital 
computer  root  locus  program  is  developed  to  determine  feedback 
coefficients  in  a  manner  which  minimizes  the  possibility  of  oscilla- 
tory response  in  the  presence  of  state  sensor  errors.  The  effect 
of  noise  on  the  state  signal  is  investigated  and  a  sample  data  fil- 
tering technique  developed.  A  quasi-optimum  time  technique  utilizing 
second  order  switching  logic  for  initial  control  effort  and  linear 
feedback  for  terminal  control  is  developed. 
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1 .0  Introduction 

The  purpose  of  this  study  was  to  investigate  means  of  controlling 
fourth  order  systems  to  obtain  various  types  of  responses.  The  initial 
portion  of  the  study  deals  with  the  use  of  linear  state  signal  feedback 
to  obtain  exponential  settling.  Use  of  the  uncoupled  form  is  investi- 
gated for  use  in  the  feedback  solution  analysis.  The  effect  of  noise 
on  the  state  signals  for  the  linear  feedback  solution  was  investigated 
and  a  sample  data  filtering  technique  developed.  A  root  locus  program 
was  developed  to  determine  feedback  coefficients  in  a  manner  which  will 
minimize  the  possibility  of  oscillatory  response  in  the  presence  of 
state  sensor  errors. 

The  remainder  of  the  study  deals  with  a  qua si -optimum  time 
solution  which  incorporates  switching  logic.  Full  control  effort  is 
used  Initially  and  switching  takes  place  at  predetermined  levels  of 
state  variable  combinations.  After  switching^  control  is  allowed  to 
decay  for  a  linear  termination  of  the  solution.  The  objective  here 
is  to  obtain  a  near-  optimum  time  response  with  no  possibility  of  con- 
trol chatter. 

Section  2.0  of  this  report  contains  the  linear  feedback  portion 
of  the  investigation.   Section  3.0  contains  the  quasi-optimum  investi- 
gation. All  synthesis  for  the  study  was  conducted  on  the  school's 
Control  Data  Corporation  160U  digital  computer  using  fortran.  All 
programming  will  be  referenced  in  the  text  and  shown  in  the  appendices. 


2.0  Linear  Feedback  for  Exponential  Settling 

The  response  of  a  system  to  a  set  of  initial  conditions  may  be  con^ 
trolled  by  feedback  of  the  system  state  variables.  That  is,  the  forcing 
function  of  the  system  Is  made  up  of  predetermined  amounts  of  each  of 
the  system's  state  variables.  Thus,  the  claaracteristic  equation  of  the 
controlled  system  may  be  adjusted  to  obtain  the  desired  response. 

If  an  exponential  response  of  the  system's  position  state  variable 
Is  desired,  the  characteristic  equation  will  have  only  negative  real 
roots,  with  a  dominant  root  that  causes  the  desired  exponential  path 
after  the  decay  of  the  non- dominant  roots.  The  proper  feedback  co=> 
efficients  of  each  of  the  state  variables  may  be  calculated  to  give  the 
desired  characteristic  equation. 
Examp le 

An  aircraft  landing  flare  is  representative  of  a  class  of  auto- 
matic control  problems  in  which  a  system  has  Initial  conditions  of  each 
state  variable  and  it  is  desired  that  the  position  state  variable 
settle  to  zero  in  an  exponential  manner.   In  order  to  design  a  suit- 
able control  system,  a  mathematical  description  of  the  aircraft  longi- 
tudinal motion  is  required.  Assuming  constant  air  speed  and  a  shallow 

glide  angle  leads  to  the  short  period  equation  of  longitudinal  motion 

|i  2] 

L '  J  .   These  are  written  in  terms  of  the  following  transfer  function 

relating  elevator  position,  d  (radians),  to  aircraft  pitch  angle, 

p  (radians); 

/  X       K(Ts  +   1) d(s)  ,,- 

p(s)  = ^ ^ ^—  ^  ^  (1^ 

s(s  +  2as  +  w  ) 


where     K  =  short  period  gain 

w  =  short  period  resonant  frequency 
a  =  short  period  damping  factor 
T  =  path  time  constant 
In  order  to  complete  the  mathematical  description  of  the  aircraft, 
the  transfer  function  relating  altitude,  h  (feet),  and  pitch  angle, 
p  (radians),  in  terms  of  velocity  V  (feet  per  second)  and  path  time 
constant  T  is 

V        --  (2) 


h(s)  = 


P(s) 


s(Ts  +  1) 

Combining  (l)  and  (2)  results  in  a  transfer  function  relating 
altitude  and  elevator  position: 


h(s)  = 


KV 


2/  2   ^      2x 
s  (s  +  2as  +  w  ) 


d(s) 


(3) 


Letting  the  system  forcing  function   u   be  equal  to  KVd,  z   signal 
flow  diagram  of  the  system  is  shown  in  Figure  1. 


Figure  1.   Block  Diagram  of  Aircraft  Longitudinal  Motion  with 
State  Variable  Feedback 


The  characteristic  equation  for  the  system  of  Figure  1  is 


s'  +  (b,  +  2a)s^  +  (b  +  w  )s  +  l^s  +  b^  =.  0  {h) 


It  is  now  possible  to  solve  for  the  feedback  coefficients 
b^ ,  b  ,  b  ,  and  b,  such  that  the  characteristic  equation  will  have 
the  desired  roots.   Suppose,  for  a  particular  aircraft 
a  =  0.5  and  w  =  1.0,  and  the  desired  closed  loop  characteristic  equa- 
tion roots  are  s  =-0.l8,  -1.0,  -1,0,  and  -5.O;   giving  the  charac- 
teristic equation 

s  +  7.l8s^  +  12.26s^  +  6.98s  +  0.9  ^  0.  (5) 

Equating  (4)  and  (5)  will  then  give  the  desired  feedback  co- 
efficients 

^1  ^   °-^'  \   "  ^'^^'   ^3  '   11«26,  b^  -  6.18.  (6) 

This  will  give  a  time  response  solution  for  the  flare^  after  allowing 
a  short  decay  time  (t^  :^  6  seconds)  for  the  non-dominant  roots^,  of 
approximate ly 

h(t)  =  h(t^)  e-0-^^*^  (7) 

2.1  Use  of  the  Uncoupled  Form 

In  order  to  more  closely  examine  sources  of  instability  in  the 
above  type  of  problem,  and  to  determine  the  proper  variable  for  a  root 
locus  study,  it  is  useful  to  solve  for  the  system' s  uncoupled  form   fs  j . 
Rewriting  (3)  in  the  time  domain  gives 

h(t)  +  2  ah(t)  +  w^  h(t)  =  KVd(t)  (8) 


In  order  to  rewrite  (8)  in  matrix  form  let 


x^(t)  =  h(t) 
Xp(t)  =  h(t) 

•  Co 

x^(t)  =  h(t) 

x^(t)  =  h(t) 
u(t)  =  KVd(t) 
This  leads  to  the  following  matrix  form  for  the  system; 


(9) 


0 

1 

0 

0 

0 

0 
0 

0 
0 

1 

°2 
-w 

0 

1 

X   + 

0 
0 

0 

0 

-2a 

1 

X  = 


which  is  defined  as 

X  =  Fx  +  Du 
The  system  can  now  be  transformed  into  the  uncoupled  form  by 
defining  a  new  variable  y  as 

yj^  =  u(s)/s 

y  =  u(s)/(s  +  2as  +  w  ) 
y,  =  u(s)s/(s  +  2as  +  w  ) 
After  expanding  x  by  partial  fraction  expansion  it  is  noted  that 


(10) 


11 


1/w^ 
0 
0 
0 


which   is   defined  as 
X  =  G^ 


.2a/wg 

1/w 

0 

0 


(Ua^ 


w   )/w 
-2a/w'' 

2a/w^ 
-1/w 

1 

0 

0 

1 
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{1\) 


Also,  by  solving  for  the  inverse  of  G,  the  expression 


2  =  G  X 
can  be  written  as 


(15) 


1  = 


w  2a  1     0 

0  w  2a    1 

0  0  10 

0  0  0     1 


(16) 


In  order  to  draw  a  signal  flow  diagram  of  the  uncoupled  system^ 
(li+)  is  substituted  into  (11)  giving 

G^  =  FG^  +  Du  (17) 

1 


and  multiplying  by  G 

^  =  G"  FG2  +  G"  Du 
Solving  (l8)  gives  the  uncoupled  form 


18) 


1  = 


0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

°2 
-w 

1 

0 

0 

-2a 

1  + 


which  is  drawn  as  shown  in  Figure  2.   Lines  are  added  to  denote 
subsequent  addition  of  state  variable  feedback  loops. 


n^. 


(^ 

/ 

^ 

^2 

/ 

1 

s 

1 
s 

^1 

_  / 

^    u 

Figure  2.   Block  Diagram  of  Uncoupled  Form  of  Aircraft 

Longitudinal  Motion  with  State  Variable  Feedback 

The  characteristic  equation  of  the  uncoupled  system  of  Figure  2  is 

s  +  (2a+n  +n|  )s^+  (w  +n.+2an  +n  )s  (20) 

P  2 

+  (2an.+w  n  )s  +  n.w  =  0. 

Assuming  the  same  aircraft  as  before  with  a  =  0.5  ^i^d  w  -  1 .0  gives 

s  +  (1+n  +n,  )s  +  (l+n.+n  +n  )s  +  (n.+n  )s  +  n^  =  0         (21) 

Suppose,  prior  to  beginning  the  flare,  the  aircraft  is  descend- 
ing a  -h     glideslope  with  an  airspeed  of  170  knots  and  that  the 
flare  will  begin  at  100  feet.   Nominal  initial  vertical  velocity  is 
approximately  -20  feet  per  second.   In  order  for  the  glide  angle  to 


be  approximately  tangent  to  the  flare  path^  a  dominant  characteristic 
equation  root  of  -0.2  is  selected.   That  is^  the  slowest  phase  plane 
eigenvector  is  placed  through  the  nominal  initial  condition  in  the 
x,x  phase  plane.   Now  the  feedback  coefficients  for  the  control  system 
can  be  calculated.   Suppose  characteristic  equation  roots  are  selected 
at  s  =  -0.2,  -1.0,  -1.0,  and  -5.0.   This  gives  the  following  desired 
characteristic  equation; 


s^  +  7.2s^  +  12. Us^  +  7»2s  +1=0 


Equating  (21)  and  (22)  gives  the  following  solution  for  the  U3i= 
coupled  state  variable  feedback  coefficients 

n,  =  1.0,  n  =  6.2,  n  =  i<-.2,  n,  =  0 

which  can  be  written  in  matrix  form  as 

N  =    1.0   6.2   k.2        0 
The  solution  for  the  original  state  variable  feedback  coeffi- 
cients can  be  calculated  by  noting 

B  =  NG" 
which  yields 


B  = 


1.0   7-2  n.h       6.2 


2 .2  Use  of  a  Root  Locus  Study 

In  order  to  become  familiar  with  the  behavior  of  the  roots  of 
the  characteristic  equation  (21),  a  fortran  root  locus  program  was 
written.   Since  the  uncoupled  characteristic  equation  contained  n 
as  an  element  of  each  of  the  internal  coefficients,  n  was  used  as 
the  variable  for  the  root  locus.  As  shown  in  Appendix  I^  it  was 


:22) 


(23) 


:2k) 


(25) 


(26) 
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found  that  {2k)   had  all  negative  real  roots  for  the  dominant  root  vary- 
ing from  -0o2  to  -O.366.  Thus  the  previous  solution  is  on  the  boundary 
where  a  slight  error  in  a  feedback  setting  could  cause  two  of  the 
roots  to  leave  the  real  axis  and  result  in  a  oscillatory  component 
in  the  system  response.   For  this  reason,  it  is  better  to  solve  for 
a  dominant  root  of  about  -0.l8  and  then  adjust  n^  to  move  the  dom- 
inant root  back  to  -0.2.   Doing  this^,  as  shown  in  Appendix  I  gives 

N  =   [0.9   5.71   ^.28   0.1  ]  (27) 

B=  [0.9   6.61  10.89   5.81]  (28) 


and 


These  feedback  coefficients  give  all  real  roots  for  a  dominant 
root  varying  from  -O.I8  to  -0.28^4-  when  n  is  varied  from  ^,96  to 
6.09"   Figure  7  shows  the  system  time  response  resulting  from  use  of 
the  state  variable  feedback  coefficients  of  (28). 
2.3  Consideration  of  Acceleration  Constraint 

The  proper  method  for  selecting  the  non- dominant  root  locations 
of  the  closed  loop  system  characteristic  equation  has  been  ignored 
in  the  previous  discussion.   It  will  now  be  shown  that  system  accelera= 
tion  during  exponential  settling  is  a  function  of  both  initial  con- 
ditions and  the  roots  of  the  closed  loop  system  characteristic  eq- 
uation. Thus,  if  the  most  severe  initial  conditions  that  the  system 
can  be  expected  to  be  subjected  to  are  predictable,  an  acceptable 
location  for  the  non-dominant  roots  can  be  determined  „ 

As  before,  the  dominant  root  should  be  selected  to  place  the 
slow  eigenvector  through  the  nominal  initial  condition  in  the  x  x 
phase  plane.   As  the  non-dominant  real  roots  of  the  closed  loop 
system  characteristic  equation  are  moved  to  the  left  from  the  origin 


in  the  s  plane,  the  acceleration  maximum  during  the  early  portion  of 
the  response  is  increased.   Therefore,  the  non-dominant  roots  can  be 
selected  by  considering  the  maximum  amount  of  acceleration  to  which 
the  system  should  be  subjected.   In  the  aircraft  flare  problem^  the 
proper  acceleration  constraint  would  be  determined  by  structural  con= 
siderations  as  well  as  passenger  comfort. 

In  order  to  illustrate  the  effect  of  non-dominant  root  placement, 
a  graphical  solution  of  the  acceleration  on  a  fourth  order  system 
during  exponential  response  was  made.  The  closed  loop  fourth  order 
system  was  generalized  by  lumping  the  feedback  coefficients  with  the 
plant  coefficients  for  each  state  variable.  This  resulted  in  a  gen- 
eralized fourth  order  system  as  shown  in  Figure  3  with  initial  con- 
ditions  on  the  state  variables. 


h^(0) 


^2^0) 


Figure  3.   Block  Diagram  of  General  Fourth  Order  System  with 
Initial  Conditions  on  the  State  Variables 
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Using  signal  flow  techniques^,  the  state  variable  transforms  can 
be  expressed  in  terms  of  the  initial  conditions  on  the  state  variables. 
As  is  often  the  case  for  mechanical  systems^,  the  magnitude  of  maximum 
control  effort  will  be  determined  by  the  acceleration  state  variable „ 
The  operational  expression  for  the  acceleration  state  is 

s\^(0)  +  h A0){i^+c^s^)   -  h2(0)(c2S+c^)  -  h^(0)sc^ 

h  (s)  =  j^ -^ — — (29) 

s  +  c,  s  +  c  s  +  c  s  +  c 

If  a  purely  exponential  response  on  the  position  state  variable  is 
desired,  the  roots  of  the  closed  loop  system  characteristic  equation 
must  be  real.  This  can  be  expressed  mathematically  as 


^    ^    2 
s  +c^s  +c^s  ^-c^s+Cj^  =  (s+Pj^)(s+P2)(s+p  )(s+p^) 


(30) 


^2^2 


^3^3 


2 


which  leads  to  the  following  time  expression  for  the  system  acceleratioi 

o " 
h^(t)  =  e^pj^  hj^(0)(p2P^Pi^)+h2(0)(p2P+P2Pi^+P3Pi^)+  ( 

h3(0)(p2+P^+Pi^)+h^(0)   + 

o  r 

h^(0)(p^P2p3)+h2(0)(pj^p^+p^p^+P3P^)+ 

h3(0)(p^+p^+p^)+hj^(0)  + 

h^(0)(p^P2P;^)+h2(0)(p^P2+PjP2^+P.pPi^)+ 

h]^(0)(pj^P2P^)+h2(0)(pj^P2+P^p3+P2p2)-J' 
h3(0)(p^+P2+P3)+h^(0) 


Jl) 
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e,  =  -; w^   ^  M r->  etc.  (32) 


where 

'l  "   (P2-Pi)(P3-Pi)(Pi+-Pi) 

Curves  may  then  be  plotted  for  various  initial  conditions  and 
characteristic  equation  roots.   Figure  ^  is  a  graphical  solution 
obtained  from  the  fortran  program  in  Appendix  II «   It  is  a  three 
dimensional  plot  of  system  acceleration  resulting  from  various 
initial  conditions  on  the  position  and  velocity  state  as  well  as 
various  non-dominant  root  locations.   The  dominant  root  location 
is  at  -0o2. 
2  .h     Statistical  Considerations 

This  portion  of  the  study  was  made  to  investigate  the  effect 
of  noisy  state  variable  signals  when  linear  feedback  is  used.   This 
would  be  the  actual  case  in  the  previous  aircraft  flare  example 
since  a  predictable  amount  of  error  due  to  the  aircraft  state  sensors 
could  be  expected.   Radio  altimeter  errors  would  be  caused  by  thermal 
gaussian  noise,  terrain  uneveness,  and  the  effects  of  close  proximity 
to  the  ground  during  the  flare.   The  use  of  an  altimeter  as  the 
prime  sensor  would  necessitate  differentiating  three  times  in  order 
to  obtain  four  state  variables.   Because  of  the  amount  of  noise  that 
could  be  expected  on  the  altimeter  output^  this  would  be  a  highly 
unlikely  approach.   Since  an  accelerometer  would  not  be  sensitive 
to  terrain  uneveness  and  proximity^,  a  simulated  accelerometer  output 
was  used  as  the  prime  system  sensor  for  the  computer  synthesis. 
A  simulated  altimeter  output  was  used  as  a  secondary  system  sensor . 

12 


Initial  Position  100  feet 

Non-dominant  Roots  at 

a.  -1.0 

b.  -2.0 

c.  -3.0 

d.  -4.0 

e.  -5.0 

f.  -6.0 


Initial 

Velocity 

(ft/sec) 


Acceleration 
(ft/sec/sec) 


20 


15 


10 


0  Iz^ 


Figure   4       Acceleration   of  Foui^th  Order  System  with  Dominant 
Characteristic   Equation  Root  of  0.2 
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Sensor  noise  was  simulated  using  a  gaussian  noise  generator  Uy. 
See  Appendix  III  for  the  mathematical  development  and  computer  check- 
out of  this  noise  generator. 

A  digital  computer  was  used  to  simulate  the  control  problem.  A 
fourth  order  Runge-Kutta  numerical  integration  algorithm  was  used  to 
simulate  system  motion. 

As  in  the  previous  example,  it  was  desired  that  the  system  be  con- 
trolled over  an  exponential  flare.   The  plant  of  Figure  1  with 
a  =  0.5  and  w  =  1 .0  was  used.   The  feedback  coefficient  of  equation 
(28)  were  selected  to  give  a  dominant  closed  loop  characteristic 
equation  root  at  S  =  -0.2. 

In  the  computer  synthesis,  the  simulated  accelerometer  output^ 
with  gaussian  noise,  is  used  to  calculate  the  four  system  states. 
Smoothing  is  used  to  minimize  the  effects  of  the  added  noise. 
Periodically,  the  altitude  state  is  updated  with  a  simulated  noisy 
altitude  output  which  has  also  been  smoothed.   Updating  is  accomplish- 
ed by  averaging  the  two  altitudes. 

Smoothing  of  both  the  simulated  acceleration  and  altitude  is 
accomplished  by  calculating  a  least  square  error  line  over  a  vari- 
able number  of  past  sensor  outputs.   See  Appendix  IV  for  a  mathe- 
matical development  of  the  smoothing  method. 

Figure  5  shows  the  time  response  of  the  system  in  a  three 
dimensional  graph.   This  multi-curve  graph  shows  the  effects  of 
accelerometer  noise  variation  on  the  response.  Eleven  curves  are 
shown  with  accelerometer  noise  variance  varying  from  zero  to  0,01 
feet  per  second  per  second.   Lines  joining  points  of  equal  time 

Ik 


Altitude 
(feet) 


100 


Variance  of 
AcceleroXiGter  Noise 
.001  feet/sec/sec) 


Figure  5   Time  Response  of  Fourth  Order  System  with 
Noisy  Sensors 
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(l  second^  2  seconds,  etc.)  are  shown  to  give  an  illusion  of  depth  to 
the  graph.   For  each  curve  the  altitude  noise  variance  was  0.25  feet. 
Linear  smoothing  over  three  sample  pointS;,  having  a  time  duration  of 
0.025  second  apart,  was  used  for  both  acceleration  and  altitude. 
These  curves  show  that  the  terminal  portion  of  the  response  becomes 
erratic  when  accelerometer  variance  is  0.003  feet  per  second  per 
second  or  more . 

Figure  6  is  similar  to  Figure  5  except  that  the  third  axis  shows 
the  effect  of  using  a  variable  number  of  past  sample  points  for  linear 
smoothing.   Eleven  curves  are  shown  with  the  number  of  sampling 
points  varying  from  zero  (no  smoothing) to  «-welve  ,  Acceleration 
variance  was  O.OO5  feet  per  second  per  second  and  altitude  variance 
was  0.25  feet.   It  is  seen  that  as  sampling  points  are  increased  to 
five  or  more,  the  data  becomes  sufficiently  stale  so  that  oscillation 
occurs.  As  the  number  of  past  sample  points  is  increased  even  more>, 
instability  occurs. 
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Figure  6 


Time  Response  of  Fourth  Order  System  with  Noisy 
Sensors 
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3.0  Use  of  Second  Order  Switching  Logic  for  Quasi- 
Optimal  Response 

Minimum  time  solutions  for  specific  fourth  order  systems  can  be 
found  in  the  literature   [5]  «  These  solutions  are  derived  through 
the  use  of  Pontryagin  s  maximum  principle  „ 

At  this  time,  a  general  solution  for  the  fourth  order  minimum 
time  problem  has  not  been  developed  because  of  the  high  degree  of 
complexity  involved.   This  portion  of  the  investigation  was  conducted 
to  investigate  the  application  of  the  well  developed  second  order 
minimum  time  solutions  to  the  fourth  order  problem  60 

The  time  optimal  solution  for  fourth  order  systems  requires 
continuous  maximum  control  ef fort „   For  all  but  discrete  sets  of 
initial  conditions^,  three  switching  points  where  control  effort  re- 
verses, are  required „   In  this  study,  a  quasi-optimal  solution  was 
obtained  by  using  switching  logic  for  only  one  switching  point.  The 
philosophy  is  to  use  maximum  control  effort  during  the  initial 
portion  of  the  solution  and  linear  control  for  exponential  settling 
during  the  terminal  portion  „   Switching  logic  is  used  to  terminate 
the  maximum  control  effort  at  the  proper  time.   This  switch  must 
take  place  sufficiently  early  to  prevent  excessive  overshoot^  yet 
late  enough  to  result  in  a  reasonably  fast  response.   Control  effort 
is  allowed  to  decay  during  the  terminal  portion  of  the  solution.   This 
linear  portion  then  replaces  the  last  two  switches  of  the  minimum  time 
solution  and  therefor  prevents  chatter  mode. 
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In  section  2.1  it  was  shown  that  the  fourth  order  system  could  be 
represented  in  the  uncoupled  form  as  two  second  order  systems .   E7 
doing  this^  the  well  developed  second  order  switching  lines  of  Titus 
and  Demetry  can  be  utilized  to  give  an  approximate  type  of  switchirug 
logic  for  the  fourth  order  system. 

Since  two  second  order  systems  result  from  the  uncoupling  of  the 
fourth  order  system,  some  sort  of  priority  must  be  used  to  determine 
which  uncoupled  state  pair  should  be  applied  to  the  second  order  switch- 
ing criteria.   Various  priority  schemes  were  used  includiimg  selection, 
of  the  second  order  uncoupled  state  pair  having  higher  velocity^  higher 
energy,  and  variations  of  the  latter.   Also  a  switching  line  which 
averaged  the  states  of  both  the  uncoupled  state  pairs  was  tried  „ 

The  mathematical  expression  for  minimum  time  second  order  switch- 
ing lines  varies  according  to  the  type  of  plant  eigenvalues  involved 

\6j.      For  this  investigation,  the  second  order  switching  line  for 
null  roots  was  used.  This  switching  line  is 


control  =  -N  sgn 


X,  + 


h\h 


I  2N 

where  N  is  the  saturated  control  effort. 

A  digital  computer  was  used  to  investigate  this  type  of  control 
See  Appendix  V  for  the  digital  computer  program.   In  each  of  the 
following  cases,  the  plant  of  Figure  1  with  a  =  0.5  ^nd  w  =  1.0  was 
used.  The  terminal  linear  control  was  determined  by  the  feedback 
coefficients  described  in  equation  (28)  which  give  a  dominant  charac- 
teristic  equation  root  of  -0.2. 
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Figure  7  shows  the  system  time  response  when  only  saturating 
linear  control  is  used.   Initial  system  position  for  each  curve  is 
100  feet.   Initial  system  velocity  is  varied  frc»-m  30  feet  per  second 
to  -30  feet  per  second.   Control  saturation  occurs  at  plus  or  minus 
100.  Time  required  for  completion  of  90  percent  of  the  desired  travel 
varies  from  approximately  7  seconds  to  17  seconds  for  the  extreme 
initial  conditions.  These  travel  times  will  be  used  to  evaluate  the 
following  control  systems  which  use  full  control  during  the  initial 
portions  of  each  solution. 

Figure  8  shows  the  system  time  response  when  the  second  order 
switching  logic  of  equation  (33)  is  applied  to  the  average  of  the 
uncoupled  positions  and  the  average  of  the  uncoupled  velocities. 
The  initial  control  equation  therefore  is  changed  to 


control  -  "N  sgn 


13  2N 

Initial  conditions  and  control  saturation  remain  the  same ^  Time  re^ 

quired  for  completion  of  90  percent  of  the  desired  travel  varies  from 
approximately  3  seconds  to  10  seconds.  This  switching  logic  causes 
the  system  to  respond  approximately  twice  as  fast  as  it  did  when  only 
linear  control  was  used.  However,  this  control  logic  allows  overshoot 
when  the  system  has  negative  initial  velocity o   Therefore^  it  may  be 
unacceptable  for  some  applications. 

Figure  9  shows  the  system  time  response  when  the  second  order 
switching  logic  of  equation  (33)  is  applied  to  the  uncoupled  state  pair 
having  the  higher  velocity.   Initial  conditions  and  control  saturation 
level  remain  the  same.   This  shows  that  the  use  of  the  uncoupled 
velocity  to  determine  the  uncoupled  state  pair  priority  is  not  accept- 
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able.  For  initial  system  velocities  of  10  feet  per  second  or  more,  the 
response  is  identical  to  the  linear  response  of  Figure  7o  For  initial 
system  velocities  of  zero  feet  per  second  or  less^  excessive  overshoot 
occurs . 

Figure  10  shows  the  system  time  response  when  the  second  order 
switching  logic  of  equation  (33)  is  applied  to  the  uncoupled  state 
pair  having  the  higher  energy.  To  do  this,  each  of  the  uncoupled 
velocity  states  was  normalized  by  dividing  by  system  natural  frequency, 
w.   Then  the  uncoupled  state  pair,  whose  states  were  farther  from  the 
phase  plane  origin,  were  assumed  to  have  the  higher  energy.   Initial 
conditions  and  control  saturation  level  remain  unchanged.  A  fairly 
consistent  overshoot  resulted  from  all  the  initial  conditions  con= 
sidered.  Therefore,  this  control  logic  apparently  has  some  potential. 
However,  the  logic  must  be  altered  to  cause  the  switching  point  to 
occur  sooner . 

Figure  11  is  similar  to  that  of  Figure  10  except  that  the  control 
logic  was  changed  to  cause  the  switching  point  to  occur  earlier.   The 
initial  control  equation  was  changed  to 


control  =  -N  sgn 


Y,   + 


Y   Y 

2  I  2 


1         N 


(35) 


A  considerable  improvement  was  realized  by  changing  the  control  logic. 
The  overshoot  resulting  from  all  initial  conditions  is  now  much 
smaller.  Time  required  for  completion  of  90  percent  of  the  desired 
travel  varies  from  approximately  3  seconds  to  5  seconds  for  the  extreme 
initial  conditions. 
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Figure  12  is  similar  to  Figure  11  except  the  initial  control  logic 
was  again  changed  to  cause  the  switching  point  to  occur  earlier.  The 
initial  control  equation  was  changed  to 

Y   Y  I 
2   2  I 
control  =  -N  sgn   Yj^  +   q^^  ^^ —  .  (36) 

Initial  system  positions  are  100  feet,  70  feet^  and  i+0  feet.   Initial 
system  velocities  are  30>  0,  and  -30  feet  per  second.  Moderate  over- 
shoot occurs  for  small  initial  position  combined  with  large  negative 
velocity.   Otherwise,  the  system  is  not  extremely  sensitive  to  small 
Initial  positions.   This  control  system  would  therefore  be  acceptable 
for  some  applications  where  a  small  amount  of  overshoot  could  be 
accepted  under  extreme  initial  conditions. 
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U.O  Conclusions 

By  selecting  the  proper  closed  loop  system  characteristic  equa- 
tion roots,  a  fourth  order  system  may  be  controlled  to  give  a  desired 
exponential  response  within  initial  condition  limitations.   The 
dominant  characteristic  equation  root  will  completely  define  the 
system  output  trajectory  after  decay  of  the  components  associated 
with  the  non-dominant  roots.   Maximum  system  acceleration  is  a 
function  of  initial  conditions  and  characteristic  equation  roots „ 
A  root  locus  study  may  be  used  to  select  characteristic  equation 
roots  which  will  minimize  the  possibility  of  oscillatory  response 
components  in  the  presence  of  state  sensor  errors. 

Within  initial  condition  limitations,  quasi-optimum  time  con- 
trol of  a  fourth  order  system  is  possible  using  second  order  switch- 
ing logic  combined  with  terminal  linear  control o 

It  is  suggested  that  further  studies  be  conducted  to  investigate 
the  use  of  other  types  of  switching  logic  for  optimum  control  of 
fourth  order  system. 


29 


BIBLIOGRAPHY 

1.  Perkins^  C=  D,  and  Hage,  R,  E,  Airplane  Performance  Stability 
and  Control.   John  Wiley  and  Sons,  19^9. 

2.  Ellert,  F.  J.  and  Merriam,  C.  W.   Synthesis  of  Feedback  Controls 
Using  Optimization  Theory  -  An  Example.   IEEE  Transactions 

on  Automatic  Control,  v  3,  April  I963:   89- 103 „ 

3.  Electronics  Laboratories,  Stanford  University.   The  Analysis 
and  Synthesis  of  Nonlinear  Continuous  and  Sampled  Data  Systems 
Involving  Saturation,  by  F.  Kurzweil,  Jr ,   November^  1959' 
Technical  report  No.  2101-1. 

k.     Hamming,  R.  W.   Numerical  Methods  for  Scientists  and  Engineers, 
McGraw-Hill  Book  Company,  I962. 

5.  Stanford  University.   Optimum  and  Quasi-Optimum  Control  of 
Third  and  Fourth  Order  Systems,  by  I.  Flugge-Lotz  and  H.  A. 
Titus,  Jr.   October,  I962.   Technical  report  No.  13^. 

6.  United  States  Naval  Postgraduate  School.   Optimal  Control 
Systems  Part  II  (Synthesis),  by  Harold  Titus  and  James  Demetry. 
September,  19^3 •  Technical  report  No.  38 • 


30 


APPENDICES 


31 


APPENDIX  I 
ROOT  LOCUS  PROGRAM 


Program  RTLOCUS 


enter 


read 
data 


calculate 

polynomial 

and  call 

RTPLSUB 


sort  roots 
into  proper 
bins 
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..JOB»HARRIS  GRAPH  BOX  228 
PROGRAM  RTLOCUS 

1  visoKcoNvlso)^^^''**'^^^^^^^'^^''^^^^  **•**' '^^'^'^^^'^'^°''>^^^°^» 

^A^^i^^lLCU^^l]^^^^^^^'^^^    LOOPS, 9X,25HCHARACTERISTIC    EOUMION. 
-I    lS£5^?°?LU'iM^^°°T2,  nx,5HR0nT3,  nX,5HR00TU/2X,     ^^^^  '  *""'» 

112^n^Vi  ^SbL^*^^^  ^?b'iA^^'2HSU.3X,2HS3,k.2HS2,3X.^HSl,3X, 

PRINT  5 
10  FORMAT(I5,8F5.3) 

S      ^ii'z'?)  ^^^    COEFFS  OF  POLYNOMIAL 

C      FBLP(]-4)  ARE  FEEDBACK  LOOP  COEFFS 

r      RfV^^T'*K.^'^.^/T'^?.^'^^'^^^  I^  FEEDBACK  LOOP  COEFFS  FOR  EACH  ITERAT 
r      Tii7^^.?^?^>XU7^  /t^f  ?f^'-.^ND  IMAG  PARTS  OF  UNSORTEO  R30TS 

C      THE  SORTED  ROOTS  •■  ^'^'^^^  '^^^  ^^^^    ^^^^^    °^ 

DO  20  K=l ,800 

A(1)=1.0 

A(2)=1.0+FBLP2+FBLPU 

A(3)=1.0+FRLP1+FBLP2+FBLP3 

A(U)=FBLP1+FBLP2 

A(5)=FBLP1 

CALL  RTPLSUB  (M , A, U, V, CONV) 
,^  IF  (K-n  12,12,io6 
12  TRP(K,n=U(l  ) 

TIP(K,n=V(l) 

TRP(K,2)=U(2) 

TIP(K,2)=V(2) 

TRP(K,3)=U{3) 

TIP(K,3)=V{3) 

TRP(K,U)=U(U) 

TIP(K,U)=V(U) 

GO  TO  14 
600  DO  605  N=l, U 

DO  605  J=1,U 

i?vyN?-TiP{K-i^jTn*i^i^'^^'*^^^'^^"^^^*'^"^'^'^*^v^'^»"'^iP''<-i»J>J 

605  CONTINUE      ' 
. DO  60  1=1,4 
IT=1 

DO  1000  N=1,U 
00  1000J=1,U 
BIN(IT)=DEV(N,J) 
IT=IT+1 
1000  CONTINUE 
LAG=1  ■ 
DO  1020  L=2,16 

10,0  liM^2^?i;]T^'^"-"  '020.1020. ,0,0 

BIN( 1 )=BIN(L) 
BIN(L)=TEMP 
LAG=L 
1020  CONTINUE 

30  N=/°  '20,31.32,35,3U,35,36,37,38,39,U0,U1.42,U3,UU,U5).LAG 
J=l 

DEV(l,n  =  1000. 
DEV(1,2)=1000.  • 
DEV(1,3)=1000. 
DEV( 1,U)=1000. 
0EV(2,1 )=1000. 
DEV(3,1 )=1000. 
D£V(U,1)=1000. 
GO  TO  50 

31  N=l 
J  =  2 
DEV(1,1)=:1000. 
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32 


33 


3U 


35 


36 


37 


38 


1) 
2) 
3) 

k) 


DEV( 1,2) 
DEV(1,3) 
DEV( 1,4) 
DEV(2,2) 
0£V(3,2) 
DEV(U,2) 
GO  TO  50 

J  =  3 

DEVd 

OEVd 

0£V{  1 

DEVd 

0EV(2,3) 

DEV(3,3) 

DEV(U,3) 

GO  TO  50 

N=1 

J  =  «4 

DEV(1,1) 

DEVd, 2) 

0EV(1,3) 

DEV( 1,U) 

DEV(2,U) 

DEV(3,U) 

DEV{U,U) 

GO  TO  50 

N=2 

J=l 

DEV(2,1 ) 

DEV{2,2) 

DEV(2,3)^ 

DEV{2,U) 

DEVd, 

DEV(3, 

OEV(U 

GO  TO 

N=2 

J  =  2 

DEV(2 

DEV(2,2) 

DEV(2,3) 

DEV(2,4)= 

OEV( 1,2)= 

DEV(3,2)= 

DEV(4,2)= 

GO  TO  50 

N=2 

J  =  3 

DEV(2,n  = 

OEV(2,2)= 

DEV(2,3)= 

0EV(2,U)= 

DEV( 1,3)= 

DEV(3,3)= 

DEV(U,3)= 

GO  TO  50 

N=2 

J  =  U 

DEV(2,1 )= 

OEV(2,2)= 

DEV(2,3)= 

DEV(2,»4)  = 

DEVd,U)  = 

DEV{3,i*)  = 

DEV(i+,U)  = 

GO  TO  50 

N=3 

J=l 

DEV(3,1)= 


1)^ 

,1)  = 

50 


n  = 


1000. 
1000. 
1000. 
1000. 
1000. 
1000. 


^1000. 

nooo. 

1000. 
1000. 
1000. 
1000. 
1000. 


1000. 
1000. 
1000. 
1000. 
1000. 
1000. 
1000. 


1000. 
1000. 
1000. 
1000. 
1000. 
1000. 
1000. 


1000. 
1000. 
1000. 
1000. 
1000. 
1000. 
1000. 


1000. 
1000. 
1000. 
1000. 
1000. 
1000. 
1000. 


1000. 
1000. 
1000, 
1000. 
1000. 
1000. 
1000. 


1000. 


3h 


DEV 
DEV 
DEV 
DEV 
DEV 
DEV 
GO 

39  N  =  3 
J  =  2 
DEV 
DEV 
DEV 
DEV 
DEV 
DEV 
DEV 
GO 

hO  N=3 
J  =  3 
DEV 
DEV 
DEV 
DEV 
DEV 
DEV 
DEV 
GO 
N=:3 
J=U 
DEV 
DEV 
DEV 
DEV 
DEV 
DEV 
DEV 
GO 

U2  N  =  «+ 
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(3,2)^ 
(3,5)^ 
(3,U)^ 
(1,1)  = 
(2,1)  = 
(U,l)  = 
TP  50 


(3,1)= 
(3,2)= 
(3,3)= 
(3,U)  = 
(1,2)= 
(2,2)  = 
(U,2)  = 
TO  50 


(3.1)  = 

(3.2)  = 

(3.3)  = 
(3,U)  = 
(1,3)  = 
(2,3)  = 
(i+,3)  = 
TO  50 


(3,1)= 
(3,2)  = 
(3,3)= 
(3,U)  = 
(  1,U)  = 
(2,U)= 
(U,U)  = 
TO  50 


J=l 

DEV 
DEV 
DEV 
DEV 
DEV 
DEV 
DEV 
GO 

hZ  N  =  U 
J  =  2 
DEV 
DEV 
DEV 
DEV 
DEV 
DEV 
DEV 
GO 

UU  N=U 
J  =  3 
DEV 
DEV 
DEV 
DEV 
DEV 
DEV 
DEV 
GO 

US  N=U 


(U,l) 
(U,2) 
(4,3)  = 
(U,U)  = 
(1,1) 
(2,1)  = 
(3,1)  = 
TO  50 


(U,l)  = 
(U,2)  = 
(U,3)  = 
(U,U)  = 
(1,2)  = 
(2,2)  = 
(3,2)= 
TO  50 


{U,l)  = 
{U,2)= 
{U,3)= 
(4,U)  = 
(1,3)= 
(2,3)= 
(3,3)= 
TO  50 


1000. 
1000. 
1000. 
1000. 
1000. 
1000. 


1000. 
1000. 
1000. 
1000. 
1000. 
1000. 
1000. 


1000. 
1000. 
1000. 
1000. 
1000. 
1000. 
1000. 


1000. 
1000. 
1000. 
1000. 
1000. 
1000. 
1000. 


1000. 
1000. 
1000. 
1000. 
1000. 
1000. 
1000. 


1000. 
1000. 
1000. 
1000. 
1000. 
1000. 
1000. 


1000. 
1000. 
1000. 
1000. 
1000. 
1000. 
1000. 


J  =  U 
DEV(U,1)=1000. 
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DEV(4,2)=1000. 

DeV(U,3)=1000, 

0EV(U,i+)  =  1000. 

DEV(1,U)=1000. 

DEV(2,U)=1000. 

OEV(3,U)=1000. 

GO  TO  50 
50  TRP(K,J)=U(N) 

TIP(K,J)=V(.M) 
60  CONTINUE 

14  CONTINUE 

15  FORMAT( 1X,UF5.2,3X,5F5,2,UX,8F8.3) 

PRINT    15,FBLPl,FBLP2,FBLP3,FBLPi|,A(l  ),A(2),A(3),A(»+),A(S), 
lTRP(K,n,TIR(K,  1),TRP(K,2),TIP(K,2),  TRP{K,3)  ,TIP(K,3), 
2TRP(K,U),TIP(K,4) 

FBLP1=FBLP1+DEL1 

FBLP2  =  FBLP2  +  DEL2. 

FRLP3=FBLP3+D£L3 

FBLP»+  =  FBLP»*  +  0EL4 

NUMPTS=K 

20  CONTINUE 

DO  22  L=1,M 

DO  21  K=1,NUMPTS 

TRP(K)=TRP(K,L) 

TIP(K)=TIP(K,L) 

21  CONTINUE 

CALL  GRAPH  ( NUMPTS, TRP, T IP, 8 ) 

22  CONTINUE 
END 

SUBROUTINE  RTPLSUB ( N, A, U , V, CONV ) 

DIMENSION  A(50) ,H(50),B{50) ,C(50),Dt  50),E{50),U(50),V{5D ) ,CONV( 

F=10.0 

L  =  25 

IF(N)  52,54,52 
54  PAUSE  30 
52  NP3=N+3 

100  B(2)=0.0 
B(1)=0.0 
C(2)=0.0 
C(1)=0.0 
D{2)=0.0 
E(2)=0.0 
H(2)=0.0 

DO  101  J=3,NP3 

101  H(J)=A(J-2)  • 
T=1.0 

SK=10.0»«F 

150  IF(H(NP3))  200,151,200 

151  U(NP3)=s0.0 
V(NP3)=0.0 
C0NV(NP3)=SK 
NP3=NP3-1 
IF(NP3)152,152,150 

152  STOP  152 

200  IF(NP3-3)205,51,201 

205  STOP  205 

201  PS=0.0 
QS=0.0 
PT=0.0 
QT=0.0 
S  =  0.0 
REV=K0 
SK=10.0»»F 
IF{NP3-4)206,202,203     • 

206  STOP  206 

202  R=-H(4)/H(3) 
GO  TO  500 

203  DO  207  J=3,NP3 
IF{H(J) )204,207,204 

204  S=S+LOGF(ABSF(H(J)) ) 
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207 


208 
210 
250 


251 
252 
253 


25U 
255 

256 

257 

258 
300 
350 

351 

352 
353 
356 
35U 
UOO 
UOl 
k03 
U50 

451 

U52 
U60 
»i53 


U5U 


455 

1+56 
U57 

458 
490 


491 
492 

500 
501 


CON 
FPN 
S=E 
DO 
H(J 
IF( 
Ts- 
M={ 
DO 
S  =  H 
JJ= 
H(  J 
H{  J 
IF( 
P  =  P 
Q=Q 
GO 
HH2 
IF( 
Q=l 
P=- 
GO 
0=H 
P=( 
IF 
R  =  0 
DO 
DO 
B(  J 
C(J 
IF( 
IF{ 
AVH 
IF( 
B(N 
IF( 
AVH 
IF( 
DO 
D(  J 
E(  J 
IF( 
AVH 
IF( 
CC2 
CC3 
C(N 
CCl 
S=C 
IF( 
P=P 
Q=Q 
GO 
P=P 
Q=Q 
IF( 
R=R 
GO 
R=R 
CON 
PS  = 
QS  = 
PT  = 
QT= 
IF( 
SK  = 
REV 
GO 
IF( 
R=l 


TINUE 

1=N+1 

XPF(S/ 

208J=3 

)=H(J) 

ABSF{H 

T 

NP3-4) 

251  J  = 

(J) 

NP3-J+ 

)=H( JJ 

J)=S 

QS)  25 

S 

S 

TO  300 

=H{NP3 

HH2)  2 

.0 

2.0 

TO  257 

{NP3)/ 

H(NP3- 

(NP3-5 

.0 

490  1  = 

351  J  = 

)=H( J) 

)=B(J) 

H{NP3- 

8(NP3- 

B1=ABS 

AVHBl- 

P3)=H( 

B(NP3) 

B2=ABS 

SK-AVH 

451  J  = 

)=H{J) 

)=D(J) 

D{NP3) 

D3=ABS 

SK-AVH 

=C(NP3 

=C(NP3 

P3-1 )= 

=C{NP3 

C2»CC2 

S)455, 

-2.0. 

•{Q+1 . 

TO  456 

+(B(NP 

■»-{-B{N 

E(NP3- 

-1.0 

TO  490 

-0(NP3 

TINUE 

PT 

OT 

P 

Q 

REV)49 

SK/10. 

=  -REV 

TO  250 

T)501, 

.0/R 


FPNl  ) 

,NP3 

/S 

(4)/H(3) )-ABSF(H(NP3-l)/H{MP3) ) ) 250, 252, 252 

/2  +  3 

3,M 

3 
) 

3,254,253  .      - 


-2) 
56,255,256 


HH? 

1 )-Q«H{NP3-3) )/HH2 

)258,550,258 


1,L 
3,N 
-P« 
-pe 
D) 
1)  ) 
F{H 
SK) 
NP3 
)40 
F(H 
B2) 
3,N 
+  R* 
+  R« 
)45 
F(H 
D3) 
-2) 
-3) 
-P» 
-1  ) 
-CC 
454 

0) 


P3 

B(J-1 )-Q»B{ J-2) 

C( J-1 )-Q»C( J-2) 

352,400,352 

353,400,353 

(NP3-1 )/B(NP3-l } ) 

450,354,35U 

)-Q»B(NP3-2) 

1,550,401 

(NP3)/B(NP3)  ) 

550,450,U50 

P3 

D(J-l) 

E(J-l) 

2,500,452 

(NP3)/D{NP3) ) 

500,453,453 

CC2-Q»CC3 

1«CC3 
,455 


3-1 )*CC2-B(NP3)*CC3)/S 
P3-1 )*CC1+B(NP3)*CC2)/S 
1 ) )458,457,458 


)/E(NP3-l ) 


1,492,492 
0 


502,502 
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502 


503 

550 
551 

552 

560 
553 


55U 


555 

556 
557 


561 


558 
51 


NP=NP3-3 

U(NP)=R 

V(NP)=0.0 

CONV(NP)=SK 

NP3=NP3-1 

DO  503  J=3,NP3 

H( J)=D(J) 

IF{NP3-3)300,51.300 

IF(T)551,552,552 

P=P/Q 

Q=1.0/Q 

PP2=P/2.0 

QMPSQ=Q-PP2»PP2 

IF(QMPSQ)55U,554,553 

NP=NP3-3 

U(NP)=-PP2 

U(NP-1 )=-PP2 

S=SQRTF(QMPSO) 

V{NP)=S 

V(NP-1 )=-S 

GO  TO  561 

S=SQRTF(-QMPSQ) 

NP=NP3-3 

IF(P)555,556,556 

U(NP)=-PP2+S 

GO  TO  557 

U(NP)=-PP2-S 

U(NP-1 )=Q/U(NP) 

V(NP)=0.0 

V(NP-1 )=0.0 

CONV(NP)=SK 

C0NV{NP-1 )=SK 

NP3=NP3-2 

DO  558  J=3,NP3 

H( J)=B(J) 

GO  TO  200 

RETURN 

END 

END 


k    0.00  0.01  0.00  0.00  1.00  0.00  U.20  D.OO 


l.OOE+0 
2  HARRIS  BOX  228 

1  Rl 

2  R2 

2  R3 

3  R4 


1.00E■^0 


8 


8 
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APPENDIX  II 
Graphical  Solution  of  System  Acceleration 


Program  SETROOT 


enter 


read 
initial 
conditions 


Calculate 
acceleration 
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.JOB«HARRIS  BOX  228 

PROGRAM  SETROOT 

DIMENSION  ACCl (900) , ACC2 ( 900 ) , ACC3( 9 00 ) , ACCU ( 900) , ACCTOT (900) , 
1T{900) ,X(900) ,Y(900) 

DO  200  K=l,30 
5  FORMAT  (3F10.5) 

READ    5,ALTICtSINKIC,  ACCICBERKICPl,  P2,P3,PU 

C0EF1=P1»P1*( { ALT IC«P2»P3»PU)+SINKi:» (P2«P3+P2»PU+P3»PU) +ACCIC* 
1  (P2  +  P3  +  Pii)-t-BERKIC)/(  (P2-P1  )*(P3-P1  )»  (PU-P1  )  ) 

C0EF2=P2»P2«( (ALTIC«P1»P3«PU)+SINKi:»(Pl»P3+Pl«PU+P3»PU) +ACCIC* 
1  (PUP3  +  PU)+BERKIC)/{  (Pl-P2)*{P3-P2)t  {PU-P2)  ) 

C0EF3=P3*P3«( { ALriC»Pl»P2«PU)+SINKi: » ( P 1 «P2+P1»PU+P2»P4) +ACCIC* 
1 {Pl+P2+PU)+BERKIC)/( (P1-P3)*(P2-P3)*(P4-P3) ) 

COEFU  =  PU»PU«  (  (ALTIC»Pl»P2»Pi)+SINKi:«{Pl»P2  +  Pl»P3  +  P2*P3)  -^-ACCIC* 
1  (P1+P2  +  P3)+BERKIC)/  (  (Pl-Pl<)*(P2-PU)^  (P3-PU)  ) 

T(1)=0.0 

DO    100      J=l,300 

ACCl ( J)=C0EF1«EXPF(-P1«T( J)  ) 

ACC2( J)=C0EF2»EXPF(-P2»T(J)  ) 

ACC3(  J)=C0EF3«EX-PF(-P3*T(J)  ) 

ACC4{ J)=C0EFU»EXPF(-P4»T{J)  ) 

ACCTOK J)=ACC1 { J)+ACC2{ J)+ACC3( J)+a:CU{J) 

X( J)=T(J)-SINKIC»0.5-10.0 

Y( J)=ACCT0T( J )-S INK  I C»0.5«5. 0-50.0 

T(  J+1  )=T{  J)■^0.025 
100  CONTINUE 

CALL  GRAPH  {200,X,Y,3) 
200  CONTINUE  i 

END 

END 


100.       -20. 

♦l.OOE+00  +5.00E+00 

ACC  VS  TIME 
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APPENDIX  III 

Gaussian  Noise  Generator 

1 .   Algorithm 

Given  a  uniform  distribution  p(x)  =  1^  -.5^^"«5 

E(x)  =  -a   =   0.0 

Var(x)  =  \    (x^u)  p(x)dx  =  l/l2 

Let  R  =  c  ^  X 
1 

e(r)  ^  0„0 
Var(R)  -  Var(cnx)  ^  c  Var(  ^   x) 
-en  Var(x) 


and 


2 

c 


p(r)  ^  Norimal(Oj,c  )  by  the  Central  Limit  Theorem 


2.   Computer  Program  Flow  Chart 


enter 


/sum  12  rand . 
/numbers  having 
(  unif.  dist. 
\  from  -0..5 
\  to  Oc5 

I 

multiply 
by  desired 
deviation 


k2 


:..H,.4 


\?^\  %, 


Graph  1  Distribution  of  1000  Random  Numbers  versus  Theoretical 
Distribution  with  Standard  Deviation  of  Unity 


^3 


..JOB-HARRIS  BOX  228 

PROGRAM  CHKRAND 

DIMtlMSION    RAND(5000),KMAG(300),BIN(300),XLINE(500)  ,YLINE  (500) 

VAR=4.0 

START=71.3921 

NUMPTS=5000 

NUMBINS=100 

XNUMPTS=NUMPTS  > 

C      CALL  RANDOM  NUMBER  GENERATER  AND  PL^CE  IN  ARRAY 

DO  10    I=1,NUMPTS 

CALL  RANDOM  ( RAND{ I ) , START, VAR ) 

RANTOT=RANTOT  +  RA'ND(  I) 
10  CONTINUE  ■     ■  i 

C      CALCULATE  MEAN  AND  DEVIATION  • 

DEV=0.0 

TOT=NUMPTS 

RANMEAN=RANTOT/TOT 

DO  20  1=1 ,NUMPTS 
20  D£V=DEV+(RANMEAN-RANO(I) )«(RANMEAN-^AND( I) ) 

D£V=SQRTF(0EV/T0T) 
30  F0RMAT{//,6H  MEAN=  F10.6,5X,11H  DEVIATION=  flO.6) 

PRINT  30,RANM£AN,DEV 
C      SORTS  THE  ORDERED  RANDOM  NUMBERS  INTO  240  BINS,  INCREMENTS  BIN( 
C      FOR  EACH  NUMBER  THEREIN 

XMAG( 1 )=-3.0»SQRTF( VAR) 

UMB1NS=NUMBINS 

BININC=-2.0*XMAG{1 )/UMBINS 

DO  60   L=1,NUMBINS 

XMAG(L+1 )=XMAG(L)+BININC 

DIN{L)=0.0 

DO  50   I=1,NUMPTS 

IF  (RAND( I)-XMAG(L) )   50,45,U5 
45  IF  (RANOd  )-XMAG(L+l  )  )   U7,U7,50 
U7  BIN{L)=BIN(L)-t-l  .  0/ (  BI NJ  INC»XNUMPTS  ) 
50  CONTINUE 
60  CONTINUE 

CALL  GRAPH   ( 2U0, XM AG, BIN, 8 ) 
C      COMPUTE  ACTUAL  NORMAL  CURVE 

PI=i. 141593 

XDEV=SQRTF( VAR) 

ACTOR=SQRTF(2.0»PI) 

COEF=1.0/(ACTOR»XDEV) 

XLINEd  )=-3.0«XDEV 

DO  70   L=l,300 

XLINE{L+1 )=XLINE(L)+XDEV/50.0 

YLINE(L)=CDEF»EXPF(-0.5»(XLINE{L)/X3EV)i(XLINE(L)/XD£V)) 
.  70  CONTINUE 

CALL  GRAPH  ( 300 , XL  I NE , YL INE, 8 ) 

END 

SUBROUTINE  RANDOM  ( QUANT, START , VAR ) 
C      GENERATES  RANDOM  NUMBERS  WITH  NORMA_  DISTRIBUTION 

TOTALCL=0.0 

FACTOR=12.U5321 

DO  91    J=l, 12 

XNUM={START+61.9)»FACT0R 

NUMX=XNUM 

XINT=NUMX 

START=XNUM-XINT 

T0TALCL=T0TALCL+(START-0.5) 
91  CONTINUE 

QUANT=TOTALCL»SORTF(VAR) 

RETURN 

END 

END 

0  8    8* 

2  NOISE  DISTRIBUTION 
2  HARRIS  BOX  228 
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APPENDIX  IV 
Linear  Smoothing  of  Noisy  Data 

1 .   Mathematical  Development 

n  o  o 

minimize:  /  (y     )  =  d 
£_ ,    -'error' 

i 

define  line   y  =  Ax  -  B 

=  ^  (y.  »  Ax.  ^  B)^ 


2    '" 
then         d 


1 
n 


(1)  dd^  /dA  =  ^   2(y.  -  Ax.  -  B)(-x.  )  ==  0 

(2)  dd""  /dB  =  /   2(y.  -  Ax.  -  B)(-1)   ^  0 

1    ^     ^ 

^  y^  =  ny,   ^  x^  =  nx 


(2)    B  =  y  -  Ax 

(1)   A  =_      ^     ^  ' 


n  2 

-nxx 


^   x^ 
^  1   ^ 


i^5 


Program  NXFLARE 


enter 


read 

initial 

conditions 


add 
noise 


smooth 
acceleration 


calculate 

jerk 
velocity 
altitude 
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.JOB  HARRIS  BOX  228 
PROGRAM  NXFLARE 

THIS  PROGRAM  IS  TO  INVESTIGATE  LINEAR  CONTROL  WITH  NOISY  DATA. 
INITIAL  CONDITIONS  ARE  ASSUMED  KNOW^J  .   ONLY  ALTITUDE  AND  ACCELE 
SENSORS  ARE  USED.   CHARACTERISTIC  E3UAT I  ON- ROOTS  ARE  .2,2,2,'2- 

DIMENSION  X(U),XP{i+,900)  ,TP(900)  ,XJ(  900),YJ(900)  ,XN(U,900), 
1XFR0M3  (14,900)  .XFROmI  (  900  ) ,  X  1  NOI  SE  (  9D  0  ) ,  X3N0I  SE  (  900  ) ,  XR  ( '+  )  , 
2XX(20,11 ) ,YY(20, 1 1 ) 

START=71 .3921 

5  READ  3,N»NSAMP,X1VAR,X3VAR,T0,TF,DT,  (X( J  ),J  =  1,N) 
SAMP  =  NSAMP 

C      N=ORDER,  T0=INITIAL  TIME,  TF=FINAL  TIME,  DT=TIME  INC,  X=STATES 
T  =  TO 
3  FORMAT  (2I5,2F10.5/8F10.5) 
IF{N-13)  6,30,30 

6  PRINT  200, T, (X( J),J=1,N) 
L=0 

M=M+1 
NCOUNT=UO 
11  =  1 

TP(n)=T 

00  2  J=1,N 
XR(J)=X{J) 

2  xp(j,n)=x(j) 

GO  TO  20 

309  n  =  n-»-i  ; 

TP(n)=T  ! 

DO  12  J=1,N 
12  XP( J,I1 )=X( J) 

CALL  RANDOM  ( X3N0 ISE ( I  1 ) , START, X3VA^ ) 

XN3=X(3)+X3N0ISE( II ) 

CALL  RANDOM  { XlNOI SE ( 11 ), START ,X IVA^ ) 

XN1=X(1 )+XlNOISE{ II ) 

IF  {T-SAMP»DT)   50,60,60 
50  XR( 1 )=X( 1) 

XR(2)=X(2) 

XRI3)=X(3)  . 

XR(4)=X(U) 

XN(3,I1)=XN3 

XN(1,I1)=XN1  I 

C  NEXT  TWO  LINES  SET  UP  INITIAL  VALUES  FO^  INTEGRATION  IN  SMODTH  SUBR 

XFR0M3(2,I1 )=X(2) 
XFR0M3( 1 , n )=X( 1 ) 

XN1=0.0 

XN3=0.0 

GO  TO  310 
60  XN(3,I1 )=XN3 

XN(l,n)=XNl  • 

XN1=0.0 

XN3=0.0 

CALL  SMOOTH  ( NSAMP , DT, 3,  11 , TP,XN , XF10M3 ) 

JC0UNT=JC0UNT+1 

IF{ 10-JCOUNT)  80,80,70 
70  XR( 1 )=XFR0M3( 1,11) 

XR(2)=XFR0M3(2, II ) 

XR(3)=XFR0M3(3, n ) 

XR(4)=XFR0M3(U,I1 ) 

GO  TO  310 
80  JCOUNT=0 

C  USE  SMOOTH  ONLY  EVERY  TENTH  ITERATION  FDR  XI  THEN  TAKE  AVERAGE  POSI 
C  FOR  NEXT  ITERATION 

CALL  SMODTH  ( NS AMP, DT , 1 ,  1 1 , TP, XN, XOJT } 
XFROMI (11 )=XOUT 

hi 


n  CONTINUE 

1J=1 ,U),XFR0M1 ( II ) 

XR(  1  )  =  (XFR0M3{  1,  n  )-^XFRO^n{  II  )  )/2. 

XR(2)=XFR0M3(2, I  1 ) 

XR(3)=XFR0M3(3,I1 ) 

XR(4)=XFR0M3{U, n ) 

GO  TO  310 
310  IF  (TF-T)   10,20,20 

200  FORMAT  (/3X,UHTIME,6X,UHX( 1 ) ,6X,UHX(  2 ) , 6X, UHX ( 3 ) , 6X, UHX( U) ,6X, 
17HXl.\'0ISY,3X,7HX3N0ISY,3X,7hXlFRC]M3,  3X,  7HX2FR0M3,  3X,  7HX3  FR0M3  , 
23X, 7HXUFR0M3, 3X, 7HX 1  FROM 1//12F1 0.5)  ' 
100  FORMAT  (12F10.5) 
10  DO  15  1=1,11 

YJ(  I  )=XP( 1 , I)+X3VAR»U000. 

XJ(  I)=TP( I)+X3VAR«1000. 

IF  (NCOUNT-UO)   lU, 160, 160 
160  L=L+1 

YY(L,M)=YJ( I) 

XX(L,M)=XJ( I ) 

NCOUNT=0 
lU  NC0UNT=NC0UNT+1 
15  CONTINUE 

CALL  GRAPH  (n,XJ,YJ,8) 

GO  TO  5 
30  DO  35  L=l,21 

DO  32   M=1,l 1 

YJ(M)=YY(L,M) 

XJ{M)=XX(L,M) 
32  CONTINUE 

CALL  GRAPH  (  1  1,XJ,YJ,8) 
35  CONTINUE 

STOP 
20  CALL  RKUTTA  ( N, T , XR ,X , DT  ) 

T=T+OT 

GO  TO  309 

END 

SUBROUTINE  RANDOM  ( QJANT , START, VAR ) 

GENERATES  RANDOM  NUMBERS  WITH  GAUSSIAN  DISTRIBUTION 

TOTALCL=0.0 

FACT0R=12. 45321 

DO  91  J=l ,12  .    • 

XNUM=(START+61.9)*FACT0R 

NUMX=XNUM 

XINT=NUMX 

START=XNUM-XINT 

TOTALCL=TOTALCL+(START-0.5) 
91  CONTINUE 

QUANT=TOTALCL»SQRTF(VAR) 

RETURN 

END 

SUBROUTINE  RKUTTA  ( N, T, XR, X , DT ) 

DIMENSION  X(i+),  AK(U,U),XDOT(U),XC(U)  ,C(4),XR(U) 

C(1)=0.0 

C(2)=0.5 

C(3)=0.5 

C(U)=1.0 

DO  4  1  =  1  ,U 

TC=T+C{I)»DT 

DO  2  J=l  ,N 

2  XC{  J)=XR( J)+C( I )«AK( I-l  ,  J) 
CALL  DERIV  (XCXDOD 

DO  4  J=1,N 
U  AK( I, J)=DT»XDOT( J) 
DO  3  J=l  ,N 

3  X(J)=XR(  J)-t-(  AK(  1,J)+2.»AK(2,J)+2.»A<{3,J)+AK(U,J)  )/6. 
END 

SUBROUTINE  DERIV  {X,XOOT) 
DIMENSION  X(4),XD0T(4) 
XD0T(4)=-1 .0»X(3)-1.Q»X( 4)+C0NTR0L 
XD0T(3)=X(4) 


XD0r(2)=X(3) 

XOOTd  )=X(2} 

C0NTR0L=-0.9»X( 1)  -6.61«X(2)  -10.89»X(3)  -5.81«X{U) 

END 

SUBROUTIME  SMOOTH  ( NX , OX , J , 11 , X, Y, X<) 

DIMENSION  X(900) , Y ( U, 900 ) , XX( U,900 ),  A( 3, U ), ANSWER! U ) 

00  2  J=l,3 

00  2  1=1 ,U 

A(J,I)=0. 
2  CONTINUE 

FN  =  NX 
00  5  1=1, NX 

IN=I-1  • 

IT=n-IN 
A(3,U)=A(3,U)+Y(J, IT) 

A{2,»+)=A(2,i+)+Y(J,IT)«X{IT) 

A(1,U)^A(1,U)<-Y(J»IT)»X{IT)»X(IT) 

A{2,3)=A{2,3)+X( IT) 

A(2,2)=A(2,2)+X( IT)»X( IT) 

A{ 1 ,2)=A( 1 ,2)+X(  IT)«X( IT  )«X(  IT) 

A(l ,1 )=A( 1, 1 )+X( IT)»X( IT)»X(IT)»X(IT) 
5  CONTINUE 

Ad  ,3)=A(2,2) 

A(2,1)=A(1,2) 

A(3,1)=A(2,2) 

A(3,2)=A(2,3) 

A(3,3)=FN 

NUMBER=3 

CALL  J0RDAN2  { A, NUMBER, ANSWER ) 

IF  (J-2)  8,8,10 
10  XX (3, II )=ANSWER( 1 )*X( 11 )*X( II ) +ANS WE R ( 2 ) »X ( 11 )+ANSWER{3) 

XX(U,I1 )=2.*ANSWER(  1 )»X(  II  )+ANSWER(2) 

XX(2,n  )=XX(2,I1-1  )+XX(3,n  )«DX 

XX(1,I1 )=XX{ 1 ,11-1 )+XX(2,Il )»DX 

RETURN 
8  XX= ANSWER (1 )«X( II )»X(I1 ) +ANSWER( 2 ) •< ( 11 )+ANSWER(3) 

END 

END 
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APPENDIX  V 

Use  of  Second  Order  Switching  Logic 
for  Quasi-Optimum  Control 


Program  OPTFOUR 


enter 


read 

initial 

conditions 


calculate 
uncoupled 
states 


calculate 
opt  ^' mum 
control 


yes 


-V- 


calculate 

linear 
control 


/control 

equals 
sat 


no 


no 
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..JOB  HARRIS  DOX  228 
PROGRAM  OPTFOUR 
DIMENSION  X(U),X.P{4,900)  ,TP(900)  ,XJl  900),YJ(900) 

5  READ  3,N,T0,TFtDT, (X( J), J=l ,N) 
T  =  TO 

3  FORMAT  {I5,/8F10.5) 
IF  (N-13)   6,30,30 

6  CONTINUE 
11  =  1 

TP(n)=T 
DO  2  J=l ,N 

2  XP( J,I1 )=X( J) 

GO  TO  20  •  . 

309  n  =  n  +  i 

TP(I1)=T 

DO  12  J=1,N 
12  XP( J,I1 )=X( J) 
310  IF  (TF-T)  10,20,20 
10  DO  15  1=1,11 

YJ( I )=XP( 1,1  ) 

XJ(I)=TP(I) 
15  CONTINUE 

CALL  GRAPH  ( I  1 ,XJ,YJ,0) 

GO  TO  5         .  . 

30  STOP 
20  CALL  RKUTTA  (N,T,X,DT) 

T=T+DT 

GO  TO  309 

END 

SUBROUTINE  RKUTTA  {N,T,X,DT) 

DIMENSION  X(4),AK(i+,U),X00T(U),XC(U)  ,C(4) 

C( 1 )=0.0 

C(2)=0.5 

C{3)=0.5 

C(U)=1.0 

DO  U  1  =  1  ,U 

TC=T+C(I)*DT 

DO  2  J=l ,N 

2  XC{ J)=X( J)+C( I)«AK( I-l, J) 
CALL  DERIV  (XC,XDOT,T) 

DO  h    J=1,N 
U  AK( I, J)=DT»XDOT( J) 
DO  3  J=l ,N 

3  X( J)=X{J)+( AK( 1, J )+2.»AK(2, J)+2.»AK( 3,J)+AK(U,J) ) /6. 
END 

SUBROUTINE  DERIV  (X,XDOT,T) 

DIMENSION  Xlh) ,XDGT(U) 
C      SELECTS  UNCOUPLED  STATE  PAIRS  HAVIN3  HIGHER  ENERGY 
C      ITITIALIZE  FOR  NEW  CURVE 

IF  (T)  2,2,U 
2  JFLAG=0 

KFLAG=0 
k    XD0T(U)=-1.0«X{3)-1  .0*X(  J4)+C0NTR0L 

XD0T(3)=X{U) 

XD0T(2)=X{3) 

XDOT( 1 )=X(2) 

SAT=100. 
C      IF  SWITCHING  POINT  PREVIOUSLY  REACHED  USE  LINEAR  CONTROL 

IF  (JFLAG)  6,6,32 
C       SWITCHING  POINT  NOT  PREVIOUSLY  REA:hED 
C       CALCULATE  UNCOUPLED  STATE  VARIABLES 
6  Y1  =  X(  1  )  +  X(2)^-X(3) 

Y2  =  X(2)-^X(3}+X(4) 

Y3=X(3) 

YU  =  X(i+) 

S0MEG=1. 

EN12=S0MEG»Y1«Y1+Y2*Y2 

£N3t+=SOMEG»Y3«Y3+YU»YU 

IF     (ABSF{£N12)-ABSF(EN3U))     10,15,15 
10    Y2=Y4  • 

52 


Y1=Y3 
15  CONTINUE 
C      CALCULATE  OPTIMUM  CONTROL 

C0NTR0L=-SIGNF{SAT,Y1+Y2»ABSF( Y2)/(D.8»SAT) ) 
C      SENSE  WHEN  OPTIMUM  CONTROL  CHANGLS  SIGN 

IF  (CONTROL)  50t30,:S0 
C      CONTROL  POSITIVE  SET  JFLAG  FOR  LINEAR 
30  JFLAG=1 
100  FORMAT  (lOH  SWITCH  AT , F 1 0. 5 , UHFE ET ) 
PRINT  100, X( 1 ) 
C      CALCULATE  LINEAR  CONTROL 

.  32  CONTROL  =  -0.9»X( r)-6.61«X(2)-10.89»Xl  3)-5.81«X(4) 
C      CHECK  FOR  SATURATION 

IF  (CONTROL-SAT)  U0,U0,35 
35  CONTROL=SIGNF(SAT, CONTROL) 
i*0  KFLAG=KFLAG+1 

IF  (KFLAG-1)  45, US, 50 
U5  PRINT  200, X( 1 ) 
200  FORMAT  {15H  UNSATURATED  AT , Fl 0 .5,UHr EET ) 
50  RETURN 
END 
END  • 


0.0 

20. 

0.0?5 

UO. 

-5  0. 

+5.00E+00 

+U.00E+01 

2 

2  HARRIS 
1 

BOX 

228 

ENERGY 

0.0 

20. 

0.025 

40. 

0.0 

2 

4 

0.0 

20. 

0.025 

40. 

3D. 

2 

k 

0.0 

20. 

0.025 

70. 

-5  0. 

2 

k 

0.0 

20. 

0.025 

70. 

0.0 

2 

h 

0.0 

20. 

0.025 

70. 

3D. 

2 

U 

0.0 
2 

20. 

0,025 

100. 

-30. 

0.0 

20. 

0.025 

100. 

0.0 

2 

k 

0.0 

20. 

0.025 

100. 

3D. 

3 

55 

..END 

53 


SUBROUTINE  DERIV  (X,XDOT,T) 
DIMENSION  X(U) ,XDOT(U ) 
C      SELECTS  UNCOUPLED  STATE  PAIRS  HAVINS  HIGHER  .VELOCITY 
C      ITITIALIZE  FOR  NEW  CURVE 
IF  (T)  2,2,U 
2.  JFLAG=0 
KFLAG=0 
k    XDOT(U)=-1.0»X(3)-1.0»X(U)+CONTROL 
XD0T(3)=X(U) 
XD0T{2)=X(3) 
XD0T(1)=X(2) 
SAT=100. 
C    .  IF  SWITCHING  POINT  PREVIOUSLY  REACHiD  USE  LINEAR  CONTROL 

IF  (JFLAG)  6,6,32 
C       SWITCHING  POINT  NOT  PREVIOUSLY  REAIHED 
C       CALCULATE  UNCOUPLED  STATE  VARIABLES 
6  Y1=X(1 )+X(2)+X(3) 
Y2=X(2)+X(3)+X(4) 
Y3=X(3) 
YU  =  X{i|) 

IF  {ABSF(Y2)-ABSF(YU) )  10,15,15 
10  Y2=YU 
Y1=Y3 
. 15  CONTINUE 
C      CALCULATE  OPTIMUM  CONTROL 

C0NTR0L=-SIGNF(SAT,Y1+Y2*ABSF(Y2 )/(?.« SAT) ) 
C    •   SENSE  WHEN  OPTIMUM  CONTROL  CHANGES  SIGN 

IF  (CONTROL)  50,30,30 
C      CONTROL  POSITIVE  SET  JFLAG  FOR  LINEAR 
30  JFLAG=1 
100  FORMAT  (lOH  SWITCH  AT , Fl 0.5, UHFEET ) 
PRINT  100, X{ 1 ) 
C      CALCULATE  LINEAR  CONTROL 

32  CONTROL  =  -0.9»X{ 1 ) -6 . 6 1 »X ( 2 ) -1 0.89»Xl  3)-5.81*X(U) 
C      CHECK  FOR  SATURATION 

IF  (CONTROL-SAT)  i40,U0,35 
35  CONTROL=SIGNF(SAT, CONTROL) 
kO    KFLAG=KFLAG+1 

IF  (KFLAG-1)  45,U5,50 
45  PRINT  200, X( 1  ) 
200  FORMAT  (ISH  UNSATURATED  AT , F 1 0.5,UH= EET ) 
50  RETURN 
END 
•  END 
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SUBROUTINE  DERIV  {X,xnOT,T) 
DIMENSION  X{U),XDOT(U) 
C      AVERAGES  UNCOUPLED  VARIABLES  FOR  SWITCHING  LOGIC 
C      ITITIALIZE  FOR  NEW  CURVE 
IF  (T)  2,2,1+ 
2  JFLAG=0 
KFLAG=0 
I*  XOOT(U)=-1.0»X(3)-1  .0»X(l*)+CONTROL 
XD0T(3)=X(U) 
XD0T(2)=X(3) 
XD0T{1)=X{2) 
SAT=100. 
C      IF  SWITCHING  POINT  PREVIOUSLY  REACHED  USE  LINEAR  CONTROL 

IF  (JFLAG)  6,6,32 
C       SWITCHING  POINT  NOT  PREVIOUSLY  REAIHED 
C       CALCULATE  UNCOUPLED  STATE  VARIABLES 
6- Y1=X(1 )+X(2)+X{3) 
Y2  =  X(2)+X(3)-»-X(4) 
Y3=X(3) 
YU=X{4) 
C      CALCULATE  OPTIMUM  CONTROL 

CONTROL=-SIGNF(SAT, Y1+Y3+(Y2«ABSF( Y2 ) +Y4»ABSF( YU ) )/(2.»SAT) 
C      SENSE  WHEN  OPTIMUM  CONTROL  CHANGES  SIGN 

.  IF' (CONTROL)  50,30,30 
C      WHEN  CONTROL  CHANGES  SET  JFLAG  FOR  uINEAR 
30  JFLAG=1 
100  FORMAT  (lOH  SWITCH  AT , Fl 0. 5, 4HFEET ) 
PRINT  100, X{  1  ) 
C      CALCULATE  LINEAR  CONTROL 

32  C0NTR0L  =  -0.9«X( 1  )-6.6l*X(2)-10.89»X(  3)-5.81*X(U) 
C      CHECK  FOR  SATURATION 

IF  (CONTROL-SAT)  U0,U0,35 
35  CONTROL=SIGNF(SAT, CONTROL) 
i+0  KFLAG  =  KFLAG+1 

IF  (KFLAG-1)  U5,U5,50 
kS    PRINT  200, X(  1  )  ■ 
200  FORMAT  (15H  UNSATURATED  AT , FT  0.5 ,UH= EET ) 
50  RETURN 
END 
END 
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